{$N+,E+}  {* Link with full 8087 emulator, uses 8087 if present *}
UNIT q_FFT;

{*-----------------------------------------------------------------
	 Copyright (c) Quinn-Curtis, 1992

	 Filename:  FFT.PAS
	 Revision:  8.01    LAB
	 Date:      6/15/92

	 Description: -Fourier Analysis & Digital Filtering
-----------------------------------------------------------------*}


INTERFACE
  USES q_stdhdr, CartoCmds, tit_win;

CONST
     FFT_RECTANG = 0;
     FFT_PARZEN  = 1;
     FFT_HANNING = 2;
     FFT_WELCH   = 3;
     FFT_HAMMING = 4;
     FFT_EXACTB  = 5;


PROCEDURE WindowData(var x:   VeryLongVector;
                        nd: WORD; window:   integer);

PROCEDURE WindowFFTData(var xr:   VeryLongVector;
                        var yi:    VeryLongVector;
                        nd: WORD; window:   integer);

PROCEDURE RealFFTSolve( VAR x, sinc: VeryLongVector; n: WORD; inverse: INTEGER);
PROCEDURE FFTCalc(var xr:     VeryLongVector;
                  var yi:       VeryLongVector;
                      nd: WORD);

PROCEDURE FFTInvCalc(var xr:     VeryLongVector;
                     var yi:       VeryLongVector;
                      nd:     WORD);

PROCEDURE PowerSpectrumCalc(var xr:     VeryLongVector;
                            var yi:       VeryLongVector;
                                nd:     WORD;
                                interval:      RealType);

PROCEDURE FFT2DCalc(var xr:   recMat;
                    var yi:   recMat;
                      c,r:WORD; flag:   INTEGER);
PROCEDURE FFT2DCalc_p(xr:   precMat; yi:   precMat;
                      c,r: WORD; flag:   INTEGER);

PROCEDURE Convolve(VAR hk: VeryLongvector;
               VAR xk,yk: verylongvector;
               m,n: WORD);

PROCEDURE Freq(VAR filtcoef: VeryLongvector;
               VAR a: verylongvector;
               n: WORD; k: INTEGER);

PROCEDURE FIRFreqSample(VAR filtcoef: verylongvector;
                        fp: RealType;
                        n: WORD; dc,filtType,win: INTEGER);

FUNCTION FFTMagnitude(xr, yi : VeryLongVector; n, i :Integer) : REALTYPE;
FUNCTION FFTPhase (xr, yi : VeryLongVector; i :Integer) : RealType;
FUNCTION FFTFrequency (n : Integer; samplefreq: Realtype; i: Integer): RealType;

IMPLEMENTATION



PROCEDURE fft2 (VAR x, y: VeryLongVector; n: WORD; m, ks: INTEGER);
LABEL l1, l2, l3, l4, l5, l6;
VAR
	indx, k0, k1, k2, k3, span, j, k, kb, kn, mm, mk: INTEGER;
	rad, c1, c2, c3, s1, s2, s3, ck, sk, sq: Realtype;
	x0, x1, x2, x3, y0, y1, y2, y3: Realtype;
	cc: ARRAY[0..16] OF INTEGER;
BEGIN

 	sq := 0.707106781187;
	sk := 0.382683432366;
	ck := 0.92387953251;
        for k := 0 TO 16 DO cc[k] := 0;
	cc[m] := ks;
	mm := (m div 2) * 2;
	kn := 0;
	for k := m - 1 DOWNTO 0 DO
		cc[k] := cc[k+1] div 2;
	rad := 6.28318530718 / (cc[0] * ks);
	mk := m - 5;

	REPEAT
		kb := kn;
		kn := kn + ks;

		IF (mm <> m) THEN
		BEGIN
			k2 := kn;
			k0 := cc[mm] + kb;
			REPEAT
				DEC(k2,1);
				DEC(k0,1);
				x0 := x [k2];
				y0 := y [k2];
				x[k2] := x[k0] - x0;
				x[k0] := x[k0] + x0;
				y[k2] := y[k0] - y0;
				y[k0] := y[k0] + y0;
			UNTIL (k0 <= kb);
		END;

		c1 := 1.0;
                s1 := 0.0;
		indx := 0;
		k := mm - 2;
		j := 3;

		IF (k >= 0) THEN
			goto l2
		ELSE
			IF (kn < n) THEN
				goto l3
			ELSE
				goto l4;

l1:	   	WHILE (TRUE) DO
		BEGIN
			IF (cc[j] <= indx) THEN
			BEGIN
				indx := indx - cc[j];
				DEC(j,1);
				IF (cc[j] <= indx) THEN
				BEGIN
					indx := indx - cc[j];
					DEC(j,1);
					k := k+2;
				END
			END
			ELSE
				goto l6;
		END;

l6:		indx := cc[j] + indx;
		j := 3;

l2:		span := cc[k];
		IF (indx <> 0) THEN
		BEGIN
			c2 := indx * span * rad;
			c1 := cos(c2);
			s1 := sin (c2);
l5:
			c2 := c1 * c1 - s1 * s1;
			s2 := 2.0 * s1 * c1;
			c3 := c2 * c1 - s2 * s1;
			s3 := c2 * s1 + s2* c1;
		END;
		for k0 := kb + span - 1 DOWNTO kb DO
		BEGIN
			k1 := k0 + span;
			k2 := k1 + span;
			k3 := k2 + span;
			x0 := x[k0];
			y0 := y [k0];
			IF (s1 = 0) THEN
			BEGIN
				x1 := x[k1]; y1 := y [k1];
				x2 := x[k2]; y2 := y [k2];
				x3 := x[k3]; y3 := y [k3];
			END
			ELSE
			BEGIN
				x1 := x[k1] * c1 - y [k1] * s1;
				y1 := x[k1] * s1 + y [k1] * c1;
				x2 := x[k2] * c2 - y [k2] * s2;
				y2 := x[k2] * s2 + y [k2] * c2;
				x3 := x[k3] * c3 - y [k3] * s3;
				y3 := x[k3] * s3 + y [k3] * c3;
			END;
			x [k0] := x0 + x2 + x1 + x3; y[k0] := y0 + y2 + y1 + y3;
			x [k1] := x0 + x2 - x1 - x3; y[k1] := y0 + y2 - y1 - y3;
			x [k2] := x0 - x2 - y1 + y3; y[k2] := y0 - y2 + x1 - x3;
			x [k3] := x0 - x2 + y1 - y3; y[k3] := y0 - y2 - x1 + x3;
		END;
		IF (k > 0) THEN
		BEGIN
			k := k - 2;
			goto l2;
		END;
		kb := k3 + span;
		IF (kb < kn) THEN
		BEGIN
			IF (j = 0) THEN
			BEGIN
				k := 2;
				j := mk;
				goto l1;
			END;
			DEC(j,1);
			c2 := c1;
			IF (j = 1) THEN
			BEGIN
				c1 := c1 * ck + s1 * sk;
				s1 := s1 * ck - c2 * sk;
			END
			ELSE
			BEGIN
				c1 := (c1 - s1) * sq;
				s1 := (s1 + c2) * sq;
			END;
			goto	l5;
		END;
l3:
	UNTIL (kn >= n);
l4:
END;

{************************************************************************}

PROCEDURE rfft2 (VAR x, y: VeryLongVector;  n: WORD; m, ks: INTEGER);
LABEL l1,l2,l3;
VAR
	k0, k1, k2, k3, k4, span, j, indx, k, kb, nt, kn, mk: INTEGER;
	rad, c1, c2, c3, s1, s2, s3, ck, sk, sq: Realtype;
	x0, x1, x2, x3, y0, y1, y2, y3, re, im: Realtype;
	cc: ARRAY[0..16] OF INTEGER;
BEGIN

 	sq := 0.707106781187;
	sk := 0.382683432366;
	ck := 0.92387953251;
        for k := 0 TO 16 DO cc[k] := 0;

	cc[0] := ks;
	kn := 0;
	k4 := 4 * ks;
	mk := m - 4;

	for k := 1 TO m DO
	BEGIN
		ks := ks + ks;
		cc[k] := ks;
	END;
	rad := 3.14159265359 / (cc[0] * ks);

	REPEAT
		kb := kn + 4;
		kn := kn + ks;
		IF (m <> 1) THEN
		BEGIN
			k :=  0;
                        indx := 0;
			j := mk;
			nt := 3;
			c1 := 1.0;
                        s1 := 0;

			REPEAT
               span := cc[k];
				IF (indx <> 0) THEN
				BEGIN
					c2 := indx * span * rad;
					c1 := cos(c2);
					s1 := sin (c2);
l1: 
					c2 := c1 * c1 - s1 * s1;
					s2 := 2.0 * s1 * c1;
					c3 := c2 * c1 - s2 * s1;
					s3 := c2 * s1 + s2 * c1;
				END
				ELSE
					s1 := 0;
				k3 := kb - span;
				REPEAT
		           k2 := k3 - span;
       			    k1 := k2 - span;
					k0 := k1 - span;
		
					x0 := x[k0]; y0 := y [k0];
					x1 := x[k1]; y1 := y [k1];
					x2 := x[k2]; y2 := y [k2];
					x3 := x[k3]; y3 := y [k3];

					x [k0] := x0 + x2 + x1 + x3; y[k0] := y0 + y2 + y1 + y3;
					IF (s1 = 0) THEN
					BEGIN
						x [k1] := x0 - x1 - y2 + y3;
                        y[k1] := y0 - y1 + x2 - x3;
						x [k2] := x0 + x1 - x2 - x3;
                        y[k2] := y0 + y1 - y2 - y3;
						x [k3] := x0 - x1 + y2 - y3;
                        y[k3] := y0 - y1 - x2 + x3;
					END
					ELSE
					BEGIN
						re := x0 - x1 - y2 + y3;
						im := y0 - y1 + x2 - x3;
						x [k1] := re * c1 - im * s1;
                        y[k1] := re * s1 + im * c1;
						re := x0 + x1 - x2 - x3;
						im := y0 + y1 - y2 - y3;
						x [k2] := re * c2 - im * s2;
                        y[k2] := re * s2 + im * c2;
						re := x0 - x1 + y2 - y3;
						im := y0 - y1 - x2 + x3;
						x [k3] := re * c3 - im * s3;
                        y[k3] := re * s3 + im * c3;
					END;
					INC(k3,1);
				UNTIL (k3 >= kb);

				DEC(nt,1);
				IF (nt >= 0) THEN
				BEGIN
					c2 := c1;
					IF (nt = 1) THEN
					BEGIN
						c1 := c1 * ck + s1 * sk;
						s1 := s1 * ck - c2 * sk;
					END
					ELSE
					BEGIN
						c1 := (c1 - s1) * sq;
						s1 := (s1 + c2) * sq;
					END;
					kb := kb + k4;
					IF (kb <= kn) THEN
						goto l1
					ELSE
						goto l3;
				END;
				IF (nt =  -1) THEN
				BEGIN
					k := 2;
					goto l2;
				END;
				IF (cc[j] <= indx) THEN
				BEGIN
					indx := indx - cc[j];
					DEC(j,1);
					IF (cc[j] <= indx) THEN
					BEGIN
						indx := indx - cc[j];
						DEC(j,1);
						k := k + 2;
					END
					ELSE
					BEGIN
						indx := cc[j] + indx;
						j := mk;
					END
				END
				ELSE
				BEGIN
					indx := cc[j] + indx;
					j := mk;
				END;
				IF (j < mk) THEN
					goto l2;
				k := 0;
				nt := 3;
				kb := kb + k4;
				IF (kb > kn) THEN
					goto l3;

l2:
			UNTIL (FALSE);
l3:
		END;

		k := (m div 2) * 2;
		IF (k <> m) THEN
		BEGIN
			k2 := kn;
			k0 := kn - cc[k];
                        j := k0;
			REPEAT
				DEC(k2,1);
                                DEC(k0,1);
				x0 := x[k2];  y0 := y[k2];
				x[k2] := x[k0] - x0;
				x[k0] := x[k0] + x0;
				y[k2] := y[k0] - y0;
				y[k0] := y[k0] + y0;
			UNTIL (k2 <= j);
		END
	UNTIL (kn >= n);
END;

{************************************************************************}
PROCEDURE reorder (VAR x, y: VeryLongVector; n: WORD;  m,  ks,  reel: INTEGER);
LABEL l1, l2;
VAR
	i, j, indx, k, kk, kb, k2, ku, lim, p: INTEGER;
	temp: Realtype;

	cc, lst: ARRAY[0..16] OF INTEGER;

BEGIN

        for k := 0 TO 16 DO cc[k] := 0;
        for k := 0 TO 16 DO lst[k] := 0;
	cc [m] := ks;
	for k := m DOWNTO 1 DO
		cc [k-1] := cc [k] div 2;
	p := m - 1;
        j := p;
	i := 0;
        kb := 0;
	IF (reel=1) THEN
	BEGIN
		ku := n - 2;
                k := 0;
                WHILE (k <= ku) DO
		BEGIN
			temp := x [k + 1];
			x [k + 1] := y [k];
			y [k] := temp;
                        INC(k,2);
		END
	END
	ELSE
		DEC(m,1);
	lim := (m + 2) div 2;
	IF (p > 0) THEN
	REPEAT
		ku := cc [j] + kb;
                k2 := ku;
		indx := cc [m - j];
		kk := kb + indx;

		REPEAT
			 k := kk + indx;
   	   		 REPEAT
		         temp := x [kk];
			 	  x [kk] := x [k2];
				  x [k2] := temp;
        		  temp := y [kk];
				  y [kk] := y [k2];
				  y [k2] := temp;
				  INC(kk,1);
				  INC(k2,1);
			 UNTIL (kk >= k);

			 kk := kk + indx;
			 k2 := k2 + indx;
		UNTIL (kk >= ku);

		IF (j > lim) THEN
		BEGIN
			DEC(j,1);
                        INC(i,1);
			lst [i] := j;
			goto l1;
		END;
		kb := k2;
		IF (i > 0) THEN
		BEGIN
			j := lst [i];
			DEC(i,1);
			goto l1;
		END;
		IF (kb < n) THEN
			j := p
		ELSE
			goto l2;

l1:
    UNTIL (FALSE);
l2:
END;

{************************************************************************}

PROCEDURE trans (VAR x, y: VeryLongVector; n: WORD; eval: INTEGER);
VAR
	k, nk, nn: INTEGER;
	x1, ab, ba, y1, re, im, ck, sk, dc, ds, r: Realtype;
BEGIN
	nn := n div 2;
	r := 3.14159265359 / n;
	ds := sin (r);
	r := 2.0 * sin (0.5 * r);
	r := - r * r;
	dc := -0.5 * r;
	ck := 1.0;
	sk := 0.0;

	IF (eval=0) THEN
	BEGIN
		x[n] := x[0]; y [n] := y [0];
	END;
	for k := 0 TO nn DO
	BEGIN
		nk := n - k;
		x1 := x[k] + x[nk];
		ab := x[k] - x[nk];
		ba := y[k] + y[nk];
		y1 := y[k] - y[nk];

		re := ck * ba + sk * ab;
		im := sk * ba - ck * ab;
		y [nk] := im - y1;
		y [k] := im + y1;
		x [nk] := x1 - re;
		x [k] := x1 + re;
		dc := r * ck + dc;
		ck := ck + dc;
		ds := r * sk + ds;
		sk := sk + ds;
	END
END;

{************************************************************************}


{--------------------------------------------------------------------}

FUNCTION Parzen(n,j: RealType): RealType;
BEGIN
  Parzen := 1 - ABS((j - 0.5 * (n-1))/(0.5 * (n-1)));
END;
{--------------------------------------------------------------------}

FUNCTION Hanning(n,j: RealType): RealType;
BEGIN
  Hanning := 0.5 * (1-cos(2.0*pi*j/(n-1)));
END;
{--------------------------------------------------------------------}

FUNCTION Hamming(n,j: RealType): RealType;
BEGIN
  Hamming := 0.54 - 0.46 * cos(2.0*pi*j/(n-1));
END;
{--------------------------------------------------------------------}

FUNCTION ExactBlackman(n,j: RealType): RealType;
BEGIN
  ExactBlackman := 0.42 - 0.50 *
     cos(2.0*pi*j/(n-1) + 0.08 * cos(4.0 * pi * j / (n-1)));
END;
{--------------------------------------------------------------------}

FUNCTION Welch(n,j: RealType): RealType;
BEGIN
  Welch := 1 - SQR((j - 0.5 * (n-1))/(0.5 * (n + 1)));
END;
{--------------------------------------------------------------------}


PROCEDURE WindowData(var x:   VeryLongVector;
                        nd: WORD; window:   integer);
VAR i: INTEGER;
    multiplier: RealType;
BEGIN
  FOR i := 0 TO nd - 1 DO
  BEGIN
    CASE window OF
    0: multiplier := 1;
    1: multiplier := Parzen(nd,i);
    2: multiplier := Hanning(nd,i);
    3: multiplier := Welch(nd,i);
    4: multiplier := Hamming(nd,i);
    5: multiplier := ExactBlackman(nd,i);
    ELSE multiplier := 1.0;
    END;
    x[i] := multiplier * x[i];
  END;
END;
{--------------------------------------------------------------------}


PROCEDURE WindowFFTData(var xr:   VeryLongVector;
                        var yi:    VeryLongVector;
                        nd: WORD; window:   integer);

BEGIN
  WindowData(xr, nd, window);
  WindowData(yi, nd, window);
END;
{--------------------------------------------------------------------}

PROCEDURE RealFFTSolve( VAR x, sinc: VeryLongVector; n: WORD; inverse: INTEGER);
VAR
	 m,nn, j: INTEGER;
	 pp: RealType;
BEGIN
	m := -1;
	nn := n div 2;

	{ determine power of 2 }
	WHILE (nn > 0) DO
	BEGIN
		nn := nn shr 1;
		INC(m,1);
	END;
	nn := n div 2;

	IF (inverse = 1) THEN
	BEGIN
		pp := 0.5 / nn;
		trans (x, sinc, nn, 1);
		fft2 (x, sinc, nn, m, nn);
		for j := 0 TO nn-1 DO
		BEGIN
			x [j] := x[j] * pp;
			sinc [j] := -pp * sinc[j];
		END;
		reorder (x, sinc, nn, m, nn, 1);
		for j := 0 TO nn-1 DO		{ combine }
			x [j + nn] := sinc [j];
	END
	ELSE
	BEGIN
		for j := 0 TO nn-1 DO		{ split array in 2 }
			sinc [j] := x [j + nn];

		reorder (x, sinc, nn, m, nn, 1);
		rfft2 (x, sinc, nn, m, 1);
		for j := 0 TO nn-1 DO
		BEGIN
			x [j] := x[j] * 0.5;
			sinc [j] := sinc[j] * 0.5;
		END;
		trans (x, sinc, nn, 0);
		for j := 0 TO nn-1 DO
			sinc [j] := -sinc[j];

		for j := 1 TO nn-1 DO
		BEGIN
			x [j + nn] := x [nn - j];
			sinc [j + nn] := -sinc [nn - j];
		END
	END
END;
{--------------------------------------------------------------------}


PROCEDURE FFTSolve(VAR xr, yi: VeryLongVector;  n: WORD; inverse: INTEGER);
VAR
	m, nn, j: INTEGER;
	pp, qq: RealType;
BEGIN
	m := -1;
	nn := n;

	{ determine power of 2 }
	WHILE (nn > 0) DO
	BEGIN
		nn := nn shr 1;
		INC(m,1);
	END;

	fft2 (xr, yi, n, m, n);
	reorder (xr, yi, n, m, n, 0);
	IF (inverse=1) THEN
	BEGIN
		pp := 1.0 / n;
		qq := -pp;
		for j := 0 TO n-1 DO
		BEGIN
			xr [j] := xr[j] * pp;
			yi [j] := yi[j] * qq;
		END
	END
	ELSE
	BEGIN
		for j := 0 TO n-1 DO
			yi [j] := -yi[j];
	END
END;
{--------------------------------------------------------------------}

PROCEDURE FFTCalc(var xr:     VeryLongVector;
                  var yi:       VeryLongVector;
                      nd:     WORD);
BEGIN
   FFTSolve(xr, yi, nd, 0);
END;
{--------------------------------------------------------------------}

PROCEDURE FFTInvCalc(var xr:     VeryLongVector;
                     var yi:       VeryLongVector;
                      nd:     WORD);
BEGIN
   FFTSolve(xr, yi, nd, 1);
END;
{--------------------------------------------------------------------}

PROCEDURE PowerSpectrumCalc(var xr:     VeryLongVector;
                            var yi:       VeryLongVector;
                                nd:     WORD;
                                interval:      RealType);
VAR numdatr, timespan, normal: RealType;
    i: INTEGER;
function SquareAndSum(a,b: RealType):RealType;
begin
  SquareAndSum := sqr(a) + sqr(b);
end;

BEGIN
   numdatr := nd;
   normal := 1.0/sqr(numdatr);
   FFTSolve(xr, yi, nd, 0);
   xr[0] := SquareAndSum(xr[0],yi[0]) * normal;
   FOR i := 1 TO (nd div 2) - 1 DO
     xr[i] := normal * (SquareAndSum(xr[i], yi[i]) +
                 SquareAndSum(xr[nd-i], yi[nd-i]));
   i := nd div 2;
   xr[i] := SquareAndSum(xr[i], yi[i]) * normal;
   timeSpan := nd * interval;
   FOR i := 0 TO (nd div 2) do
     yi[i] := int(i) / timeSpan;

END;
{--------------------------------------------------------------------}

PROCEDURE FFT2DRCCalc(var xr:   recMat;
                    var yi:   recMat;
                      c,r: WORD; flag:   INTEGER);
VAR a,b,i,j: INTEGER;
    yVector,xVector: ^VeryLongVector;
label 999;
BEGIN
  yVector := NIL; xVector := NIL;

  ier := _getmem('-',pointer(yVector),sizeOf(VeryLongVector));
  if ier <> 0 then goto 999;
  ier := _getmem('-',pointer(xVector),sizeOf(VeryLongVector));
  if ier <> 0 then goto 999;
(*
  NEW(yVector);
  IF (yVector = NIL) THEN
  BEGIN
    memerrfunction(313);
    exit;
  END;
  NEW(xVector);
  IF (xVector = NIL) THEN
  BEGIN
    Dispose(yVector);
    memerrfunction(313);
    exit;
  END;
*)
  FOR j := 0 to r - 1 DO
  BEGIN
    FOR a := 0  TO c -1 DO
    BEGIN
      xVector^[a] := xr[j,a];
      yVector^[a] := yi[j,a];
    END;
    IF flag = 0 THEN
      fftCalc(xVector^, yVector^, c)
    ELSE
      fftinvCalc(xVector^, yVector^, c);
    FOR a := 0  TO c -1 DO
    BEGIN
      xr[j,a] := xVector^[a];
      yi[j,a] := yVector^[a];
    END;
  END;
  FOR i := 0 TO c - 1 DO
  BEGIN
    FOR a := 0 TO r - 1 DO
    BEGIN
      xVector^[a] := xr[a,i];
      yVector^[a] := yi[a,i];
    END;
    IF flag = 0 THEN
      fftCalc(xVector^, yVector^, r)
    ELSE
      fftinvCalc(xVector^, yVector^, r);
    FOR a := 0 TO r - 1 DO
    BEGIN
      xr[a,i] := xVector^[a];
      yi[a,i] := yVector^[a];
    END;
  END;
(*  Dispose(xVector); Dispose(yVector); *)
999:
  _freemem('-',pointer(yVector),sizeOf(VeryLongVector));
  _freemem('-',pointer(xVector),sizeOf(VeryLongVector));
END;

{--------------------------------------------------------------------}

PROCEDURE FFT2DRCCalc_p(xr:   precMat;
                        yi:   precMat;
                      c,r: WORD; flag:   INTEGER);
VAR a,b,i,j: INTEGER;
    yVector,xVector: ^VeryLongVector;
label 999;
BEGIN
  yVector := NIL; xVector := NIL;

  ier := _getmem('-',pointer(yVector),sizeOf(VeryLongVector));
  if ier <> 0 then goto 999;
  ier := _getmem('-',pointer(xVector),sizeOf(VeryLongVector));
  if ier <> 0 then goto 999;

  FOR j := 0 to r - 1 DO
  BEGIN
    FOR a := 0  TO c -1 DO
    BEGIN
      xVector^[a] := xr[j]^[a];
      yVector^[a] := yi[j]^[a];
    END;
    IF flag = 0 THEN
      fftCalc(xVector^, yVector^, c)
    ELSE
      fftinvCalc(xVector^, yVector^, c);
    FOR a := 0  TO c -1 DO
    BEGIN
      xr[j]^[a] := xVector^[a];
      yi[j]^[a] := yVector^[a];
    END;
  END;
  FOR i := 0 TO c - 1 DO
  BEGIN
    FOR a := 0 TO r - 1 DO
    BEGIN
      xVector^[a] := xr[a]^[i];
      yVector^[a] := yi[a]^[i];
    END;
    IF flag = 0 THEN
      fftCalc(xVector^, yVector^, r)
    ELSE
      fftinvCalc(xVector^, yVector^, r);
    FOR a := 0 TO r - 1 DO
    BEGIN
      xr[a]^[i] := xVector^[a];
      yi[a]^[i] := yVector^[a];
    END;
  END;
999:
  _freemem('-',pointer(yVector),sizeOf(VeryLongVector));
  _freemem('-',pointer(xVector),sizeOf(VeryLongVector));
END;
{--------------------------------------------------------------------}

PROCEDURE FFT2DCRCalc(var xr:   recMat;
                    var yi:   recMat;
                      c,r: WORD; flag:   INTEGER);
VAR a,b,i,j: INTEGER;
    yVector,xVector: ^VeryLongVector;
label 999;
BEGIN
  yVector := NIL; xVector := NIL;

  ier := _getmem('-',pointer(yVector),sizeOf(VeryLongVector));
  if ier <> 0 then goto 999;
  ier := _getmem('-',pointer(xVector),sizeOf(VeryLongVector));
  if ier <> 0 then goto 999;
{
  NEW(yVector);
  IF (yVector = NIL) THEN
  BEGIN
    memerrfunction(314);
    exit;
  END;
  NEW(xVector);
  IF (xVector = NIL) THEN
  BEGIN
    Dispose(yVector);
    memerrfunction(314);
    exit;
  END;
}
  FOR j := 0 to c - 1 DO
  BEGIN
    FOR a := 0  TO r -1 DO
    BEGIN
      xVector^[a] := xr[a,j];
      yVector^[a] := yi[a,j];
    END;
    IF flag = 0 THEN
      fftCalc(xVector^, yVector^, r)
    ELSE
      fftinvCalc(xVector^, yVector^, r);
    FOR a := 0  TO r -1 DO
    BEGIN
      xr[a,j] := xVector^[a];
      yi[a,j] := yVector^[a];
    END;
  END;
  FOR i := 0 TO r - 1 DO
  BEGIN
    FOR a := 0 TO c - 1 DO
    BEGIN
      xVector^[a] := xr[i,a];
      yVector^[a] := yi[i,a];
    END;
    IF flag = 0 THEN
      fftCalc(xVector^, yVector^, c)
    ELSE
      fftinvCalc(xVector^, yVector^, c);
    FOR a := 0 TO c - 1 DO
    BEGIN
      xr[i,a] := xVector^[a];
      yi[i,a] := yVector^[a];
    END;
  END;
{  Dispose(xVector); Dispose(yVector);  }
999:
  _freemem('-',pointer(yVector),sizeOf(VeryLongVector));
  _freemem('-',pointer(xVector),sizeOf(VeryLongVector));
END;

{--------------------------------------------------------------------}

PROCEDURE FFT2DCRCalc_p (xr:   precMat; yi:   precMat;
                      c,r: WORD; flag:   INTEGER);
VAR a,b,i,j: INTEGER;
    yVector,xVector: ^VeryLongVector;
label 999;
BEGIN
  yVector := NIL; xVector := NIL;

  ier := _getmem('-',pointer(yVector),sizeOf(VeryLongVector));
  if ier <> 0 then goto 999;
  ier := _getmem('-',pointer(xVector),sizeOf(VeryLongVector));
  if ier <> 0 then goto 999;
{
  NEW(yVector);
  IF (yVector = NIL) THEN
  BEGIN
    memerrfunction(314);
    exit;
  END;
  NEW(xVector);
  IF (xVector = NIL) THEN
  BEGIN
    Dispose(yVector);
    memerrfunction(314);
    exit;
  END;
}
  FOR j := 0 to c - 1 DO
  BEGIN
    FOR a := 0  TO r -1 DO
    BEGIN
      xVector^[a] := xr[a]^[j];
      yVector^[a] := yi[a]^[j];
    END;
    IF flag = 0 THEN
      fftCalc(xVector^, yVector^, r)
    ELSE
      fftinvCalc(xVector^, yVector^, r);
    FOR a := 0  TO r -1 DO
    BEGIN
      xr[a]^[j] := xVector^[a];
      yi[a]^[j] := yVector^[a];
    END;
  END;
  FOR i := 0 TO r - 1 DO
  BEGIN
    FOR a := 0 TO c - 1 DO
    BEGIN
      xVector^[a] := xr[i]^[a];
      yVector^[a] := yi[i]^[a];
    END;
    IF flag = 0 THEN
      fftCalc(xVector^, yVector^, c)
    ELSE
      fftinvCalc(xVector^, yVector^, c);
    FOR a := 0 TO c - 1 DO
    BEGIN
      xr[i]^[a] := xVector^[a];
      yi[i]^[a] := yVector^[a];
    END;
  END;
999:
  _freemem('-',pointer(yVector),sizeOf(VeryLongVector));
  _freemem('-',pointer(xVector),sizeOf(VeryLongVector));
END;

{--------------------------------------------------------------------}

PROCEDURE FFT2DCalc(var xr:   recMat;
                    var yi:   recMat;
                      c,r: WORD; flag:   INTEGER);
VAR a,b,i,j: INTEGER;
BEGIN
    IF flag = 0 THEN
      FFT2DRCCalc(xr, yi, c,r,flag)
    ELSE
      FFT2DCRCalc(xr, yi, c,r,flag);
END;
{--------------------------------------------------------------------}

PROCEDURE FFT2DCalc_p(xr:   precMat; yi:   precMat;
                      c,r: WORD; flag:   INTEGER);
VAR a,b,i,j: INTEGER;
BEGIN
    IF flag = 0 THEN
      FFT2DRCCalc_p(xr, yi, c,r,flag)
    ELSE
      FFT2DCRCalc_p(xr, yi, c,r,flag);
END;

{--------------------------------------------------------------------}

PROCEDURE Convolve(VAR hk: VeryLongvector;
               VAR xk,yk: verylongvector;
               m,n: WORD);
VAR nn,mm: INTEGER;
    summ: RealType;
BEGIN
  for nn := 0 to (n+m-2) do
  begin
    summ := 0.0;
    for mm := 0 to m-1 do
    begin
      IF ((nn-mm) >= 0) AND ((nn-mm) <= (n-1)) THEN
        summ := summ + hk[mm] * xk[nn-mm];
    end;
    yk[nn] := summ;
  END
END;
{--------------------------------------------------------------------}


PROCEDURE Freq(VAR filtcoef: VeryLongvector;
               VAR a: verylongvector;
               n: WORD; k: INTEGER);
VAR i,j,m,n2: INTEGER;
    att,am,q: RealType;
BEGIN
  q := pi/k;
  am := (n+1)/2.0;
  m := (n+1) div 2;
  n2 := n div 2;
  for j := 1 to k + 1 do
  begin
    att := 0.0;
    if (am = m) then att := filtcoef[m-1]/2.0;
    for i := 1 to n2 do
    begin
      att := att + filtcoef[i-1] * cos(q*(am-i)*(j-1))
    END;
    a[j-1] := 2.0 * att;
  END
END;
{--------------------------------------------------------------------}

PROCEDURE FIRFreqSample(VAR filtcoef: verylongvector;
                        fp: RealType;
                        n: WORD; dc,filtType,win: INTEGER);
VAR np,i,j,m,m1,n2: INTEGER;
    lowbin, highbin, att,xt,am,q: RealType;
    a: verylongvector;
BEGIN
  m := (n+1) div 2;
  am := (n+1)/2.0;
  m1 := n div 2 + 1;
  q := 2.0*pi/n;
  n2 := n div 2;
  if (dc = 1) THEN
    np := trunc(n * fp + 0.5)
  ELSE
    np := trunc(n * fp + 1.0);
  if filtType = 1 THEN
  BEGIN
    lowbin := 0.0;
    highbin := 1.0;
  END
  ELSE
  BEGIN
    lowbin := 1.0;
    highbin := 0.0;
  END;
  for j := 0 to np do
   a[j] := lowbin;
  for j := np+1  to m1 do
    a[j] := highbin;
  if dc = 0 then
  begin
    for j := 1 to m do
    begin
      xt := a[1] / 2.0;
      for i := 2 to m do
        xt := xt + a[i]* cos(q * (am-j) * (i-1));
      filtcoef[j-1] := 2.0 * xt / n;
      filtcoef[n-j] := filtcoef[j-1];
    end;
  end
  else
  begin
    for j := 1 to m do
    begin
      xt := 0.0;
      for i := 1 to n2 do
        xt := xt + a[i]* cos(q * (am-j) * (i-0.5));
      if (am <> m) then xt := xt + a[m]*cos(pi * (am-j))/2.0;
      filtcoef[j-1] := 2.0 * xt / n;
      filtcoef[n-j] := filtcoef[j-1];
    end;
  end;
  IF win>0 THEN WindowData(filtcoef,n,win);
end;
{--------------------------------------------------------------------}


{*--------------------------------------------------------------*}
{* Calculates the FFT magnitude associated with a given harmonic index. *}
FUNCTION FFTMagnitude(xr, yi : VeryLongVector; n, i :Integer) : REALTYPE;
var	_result, nr, rr, ii :realtype;

BEGIN
	nr := n;
	IF (i > n DIV 2) THEN
		_result := 0
	ELSE
	IF i = 0 THEN
		_result := sqrt(xr[0] * xr[0] + yi[0] * yi[0]) / nr
	ELSE
	BEGIN
		rr := abs(xr[i]) + abs(xr[n-i]);
		ii := abs(yi[i]) + abs(yi[n-i]);
		_result :=  sqrt(rr * rr + ii * ii) / nr;
	END;
	FFTMagnitude := _result;
END;
{--------------------------------------------------------------------}

FUNCTION FFTComplexAngle(c1: ComplexType):RealType;
VAR temp, CAngle: RealType;
const NearZero = 1.0e-20;
BEGIN
  with c1 do
  BEGIN
    IF abs(r) < NearZero THEN
    BEGIN
      IF abs(i) < NearZero THEN CAngle := 0;
      IF (i < -(NearZero)) THEN CAngle := 3*pi/2.0;
      IF i > NearZero THEN CAngle := pi/2.0;
      FFTComplexAngle := CAngle;
      exit;
    END
    ELSE
    BEGIN
      IF abs(i) < NearZero THEN
      BEGIN
        IF (r < -(NearZero)) THEN CAngle := pi;
        IF r > NearZero THEN CAngle := 0;
        FFTComplexAngle := CAngle;
        exit;
      END
      ELSE
      BEGIN
        IF (r > NearZero) and (i > NearZero) THEN   CAngle := arctan(i/r);
        IF (r < -NearZero) and (i > NearZero) THEN  CAngle := pi - arctan(-i/r);
        IF (r < -NearZero) and (i < -NearZero) THEN CAngle := pi + arctan(i/r);
        IF (r > NearZero) and (i < -NearZero) THEN  CAngle := 2 * pi - arctan(-i/r);
      END;
    END;
  END;
  FFTComplexAngle := CAngle;
END;

{--------------------------------------------------------------------}

{*--------------------------------------------------------------*}
{* Calculates the FFT phase associated with a given harmonic index. *}
FUNCTION FFTPhase (xr, yi : VeryLongVector; i :Integer) : RealType;
var	_result : Realtype;
    cnum : ComplexType;
BEGIN
    cnum.r := xr[i];
    cnum.i := yi[i];
	_result := FFTComplexAngle(cnum);
	FFTPhase := _result;
END;
{--------------------------------------------------------------------}

{*--------------------------------------------------------------*}
{* Calculates the FFT frequence associated with a given harmonic index
and sampling frequency. *}

FUNCTION FFTFrequency (n : Integer; samplefreq: Realtype; i: Integer): RealType;
var	_result : Realtype;
BEGIN
	IF (i <= n div 2) THEN
		_result := samplefreq * i / n
	ELSE
	BEGIN
		i := n - i;
		_result := -samplefreq * i/n;
	END;
    FFTFrequency := _result;
END;



END.
{--------------------------------------------------------------------}
